Excelencia Operacional SAG

1 Introduction

library(DT)
tags <- read.csv("tags.csv", header = TRUE, sep = ",")

options(DT.options = list(pageLength = 5))

datatable(tags) %>% formatStyle('Line', backgroundColor = styleEqual(c(0, 1, 2), c('gray', 'purple', "green"))
)

Aside from parameters fixed at design (mill dimensions, installed power, and circuit type), the major variables affecting AG/SAG mill circuit performance (throughput and grind attained) include:

1.1 Feed characteristics in terms of ore hardness/competency

The effect of feed hardness is the most significant driver for AG/SAG performance: with variations in ore hardness come variations in circuit throughput. The effect of feed size is marked, with both larger and finer feed sizes having a significant effect on throughput. With SAG mills, the response is typically that for coarser ores, throughput declines, and vice versa.

1.2 Feed size distribution

• Selection of circuit configuration in terms of liner and grate selection and closing size (screen apertures or hydrocyclone operating conditions)

• Ball charge (fraction of volumetric loading and ball size)

• Mill operating conditions including mill speed (for circuits with variable speed drives), density, and total mill load

2 Loading data

2.1 Read process data

Sub indexes 1 and 2 stand for Grind circuit 1 and 2 respectively. The operating Work index is defined as: W = \(\LARGE{ \frac{Potencia} {Feed *(\sqrt{\frac{10}{P80}} - \sqrt{\frac{10}{F80}}}}\)

dt <- as_tibble(fread("../input/conc/conc1.csv")) %>% 
  cbind(as_tibble(fread("../input/conc/df_data_interpolated_P.csv")))

# dt <- Filter(function(x) !all(is.na(x)), dt)

names(dt) <- tags$name[match(names(dt), tags$TAG)]
keep.cols <- names(dt) %in% c("NaN", NA)

### Evident conditions
dt <- as_tibble(dt[! keep.cols]) %>% 
    filter(F_M2 > 0, CC_M1>0, CC_M2>0)
###

### Data for SAG 1
dt1 <- dt %>% 
  select(c(F_M1, Csol_M1, CC_M1,
           grav1_M1, grav2_M1, peb_M1,
           hard, F80_M1, fine_1, P80_1,
           V_M1, P_M1,
           state_M1, date)) %>% 
    mutate(RPG_1 = ifelse(is.na(grav2_M1), 
                          peb_M1 + grav1_M1,
                          peb_M1 * (grav2_M1 >= grav1_M1) + 
                            (peb_M1 + grav1_M1 - grav2_M1) * 
                            (grav2_M1 < grav1_M1)),
           FF_1 = ifelse(is.na(RPG_1), F_M1, F_M1 - RPG_1), 
           F80_M1 = F80_M1 * 25.4,
           W_1 = P_M1 / (F_M1 * 
                           (10/sqrt(P80_1) - 
                            10/sqrt(F80_M1 * 1e3))
                         ),
           Ecs_1 = P_M1 / F_M1
           ) %>% 
  # filter(state_M1 == 1) %>%
  select(-c(grav1_M1, grav2_M1, peb_M1, 
            state_M1))
###

### Data for SAG 2
dt2 <- dt %>% 
  select(c(F_M2, Csol_M2, CC_M2,
           grav2_M1, peb_M2,
           hard, F80_M2, fine_2, P80_2,
           V_M2, P_M2,
           state_M2, date)) %>% 
  mutate(RPG_2 = peb_M2 + grav2_M1,
         FF_2 = F_M2 - RPG_2,
         F80_M2 = F80_M2 * 25.4,
         W_2 = P_M2 / (F_M2 * 
                          (10/sqrt(P80_2) - 
                           10/sqrt(F80_M2 * 1e3))
                       ),
         Ecs_2 = P_M2 / F_M2
         ) %>%
  # filter(state_M2 == 1) %>% 
  select(-c(grav2_M1, peb_M2, 
            state_M2))
###

### Resample time
require(reticulate)
pd <- import("pandas")

df1 <- r_to_py(dt1)
df2 <- r_to_py(dt2)

df1 = df1$set_index(pd$DatetimeIndex(df1['date']))
df2 = df2$set_index(pd$DatetimeIndex(df2['date']))

df1 = df1$resample("1H", closed = "left", label = "left", base = 8)$agg("median")
df2 = df2$resample("1H", closed = "left", label = "left", base = 8)$agg("median")

df1 = df1$dropna(how = "all")
df1 = df1$reset_index()
df2 = df2$dropna(how = "all")
df2 = df2$reset_index()

dt1 <- py_to_r(df1) 
dt2 <- py_to_r(df2)

2.2 First and last rows

2.2.1 Line 1

headTail(dt1)

##                      date    F_M1 Csol_M1   CC_M1  hard F80_M1 fine_1
## 1     2018-01-01 00:00:00 2947.63   72.59 4199.25 58.28     80  39.15
## 2     2018-01-01 01:00:00 2892.32    72.7 4197.64 58.28  72.68  42.16
## 3     2018-01-01 02:00:00 2605.67   71.13 4189.02 58.28   86.6  37.19
## 4     2018-01-01 03:00:00 2840.89   72.06 4190.68 58.28  88.65  37.52
## ...                  <NA>     ...     ...     ...   ...    ...    ...
## 14558 2019-08-31 21:00:00     0.4     NaN 1418.59  6.08    NaN  35.41
## 14559 2019-08-31 22:00:00    0.31     NaN 1403.71  6.08    NaN  33.79
## 14560 2019-08-31 23:00:00    0.54     NaN 1388.05  21.4    NaN  33.23
## 14561 2019-09-01 00:00:00    0.15     NaN 1381.33  21.4    NaN  32.37
##        P80_1 V_M1     P_M1  RPG_1    FF_1  W_1 Ecs_1
## 1     208.01 9.89  14198.4 451.48 2496.15 7.33  4.83
## 2     208.01  9.9  14329.6 482.91 2457.71 7.58  4.99
## 3     247.67  9.9 14296.54 505.06 2104.06  9.2  5.52
## 4     247.67 9.89 14393.73 100.23 2664.65  8.4  5.06
## ...      ...  ...      ...    ...     ...  ...   ...
## 14558 244.48    0        0    NaN     0.4  NaN     0
## 14559 244.48    0        0    NaN    0.31  NaN     0
## 14560 244.48    0        0    NaN    0.54  NaN     0
## 14561 244.48    0        0    NaN    0.15  NaN     0

2.2.2 Line 2

headTail(dt2)

##                      date    F_M2 Csol_M2   CC_M2  hard F80_M2 fine_2
## 1     2018-01-01 00:00:00    0.35     NaN 3540.43 58.28      0     52
## 2     2018-01-01 01:00:00 1869.03   74.43 3641.79 58.28  62.47  48.29
## 3     2018-01-01 02:00:00 1903.94   74.34 3648.57 58.28  68.13  45.93
## 4     2018-01-01 03:00:00 1907.22   74.01 3677.86 58.28  88.42   37.9
## ...                  <NA>     ...     ...     ...   ...    ...    ...
## 14558 2019-08-31 21:00:00  2712.8   75.48 3867.09  6.08 102.49  32.38
## 14559 2019-08-31 22:00:00 2291.78   75.47 3932.38  6.08 111.62   30.4
## 14560 2019-08-31 23:00:00 2702.57   75.33 3872.22  21.4 103.93  30.23
## 14561 2019-09-01 00:00:00 1914.65   71.37  3826.3  21.4  74.24  40.35
##        P80_2 V_M2     P_M2  RPG_2    FF_2   W_2 Ecs_2
## 1     211.87    0        0    NaN     NaN     0     0
## 2     211.87 6.95 10724.55  96.36  1854.8   8.4  5.44
## 3     233.82 7.56 11650.05  94.68 1761.07 10.72  6.63
## 4     233.82 8.94 14335.14 589.73 1360.19 12.24  7.61
## ...      ...  ...      ...    ...     ...   ...   ...
## 14558 236.05 9.67 14137.28 808.26 1960.44  8.43  5.22
## 14559 236.05  9.9 14303.72 713.36 1591.12 10.08  6.27
## 14560 216.72 9.89  14082.1 755.67 1591.23   8.1  5.25
## 14561 216.72  9.9 12909.11 742.33 1172.32 10.49  6.74

2.3 Dimensions

Those are the dimensions of the comminition data set:

cat('There are', c(ncol(dt1),nrow(dt1)), 'columns and rows for Line 1')

## There are 14 14561 columns and rows for Line 1

# cat('There are', c(ncol(dt2),nrow(dt2)), 'columns and rows for Line 2')

2.4 Data profile

## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning
## Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning
## -Inf

2.5 Near zero variance variables

Variables with zero variance are mostly constant across the data set, hence might provide little information and potentially cause overfitting. The table is generated on the training data with help of the caret package.

3 Missing values

3.1 Frequency

3.1.1 Line 1

gg_miss_var(dt1)

3.1.2 Line 2

gg_miss_var(dt2)

3.2 Missing data representated geometrically

The variables mentiones above are represented graphycally to search for any pattern. In this graph, the data which is missing, is given a value of 10% less than the minimum value of the available data. Then, it’s represented in a dispersion graph were the red color represents the missing data.

p1 <- dt1 %>%
  ggplot(aes(x = date, y = P80_1)) +
  geom_miss_point()

p2 <- dt2 %>%
  ggplot(aes(x = date, y = P80_2)) +
  geom_miss_point()

p3 <- dt1 %>%
  ggplot(aes(x = date, y = Csol_M1)) +
  geom_miss_point()

p4 <- dt2 %>%
  ggplot(aes(x = date, y = Csol_M2)) +
  geom_miss_point()

grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)

Fig. 1

4 Correlations

4.1 Grinding circuit Line 1

4.1.1 Histogram Line 1

dt1 %>% 
  select(-c(date)) %>%
  gather(Attributes, value) %>% 
  ggplot(aes(x=value, fill=Attributes)) +
    geom_histogram(bins = 20, colour="black", show.legend=FALSE) +
  facet_wrap(~Attributes, scales="free_x") +
  labs(x="Values", y="Frequency",
       title="Line 1 Variables - Histograms") +
  theme_bw()

## Warning: Removed 5972 rows containing non-finite values (stat_bin).

4.1.2 Pearson correlation test Line 1

4.2 Grinding circuit Line 2

4.2.1 Histogram Line 1

dt2 %>% 
  select(-c(date)) %>%
  gather(Attributes, value) %>%
  ggplot(aes(x=value, fill=Attributes)) +
  geom_histogram(bins = 10, colour="black", show.legend=FALSE) +
  facet_wrap(~Attributes, scales="free_x") +
  labs(x="Values", y="Frequency",
       title="Line 2 Variables - Histograms") +
  theme_bw() 

## Warning: Removed 10794 rows containing non-finite values (stat_bin).

4.2.2 Pearson correlation test Line 2

5 Clustering mineral hardness

5.1 Clustering

5.1.1 K-nearest neighbours

We select the variables related to mineral hardness.

# cluster_vars <- c("hard", "F80_M1", "fine_1", "P80_1") 
cluster_vars <- c("hard", "F80_M1")

dt1_knn <- dt1 %>%
  # na.omit() %>%
  drop_na(cluster_vars) %>%
  select(cluster_vars) %>%
  # select(-c(date)) %>%
  scale() %>%
  as.data.frame()

bss <- numeric()
wss <- numeric()

# Run the algorithm for different values of k
set.seed(1234)
k = 10 # Number of Clusters
for(i in 1:k){

  # For each k, calculate betweenss and tot.withinss
  bss[i] <- kmeans(dt1_knn, centers=i)$betweenss
  wss[i] <- kmeans(dt1_knn, centers=i)$tot.withinss

}

# Between-cluster sum of squares vs Choice of k
p3 <- qplot(1:k, bss, geom=c("point", "line"), 
            xlab="Number of clusters", ylab="Between-cluster sum of squares") +
  scale_x_continuous(breaks=seq(0, k, 1)) +
  theme_bw()

# Total within-cluster sum of squares vs Choice of k
p4 <- qplot(1:k, wss, geom=c("point", "line"),
            xlab="Number of clusters", ylab="Total within-cluster sum of squares") +
  scale_x_continuous(breaks=seq(0, k, 1)) +
  theme_bw()

# Subplot
grid.arrange(p3, p4, ncol=2)

5.1.2 Clustering Pair plot

isnt_out_z <- function(x, thres = 3, na.rm = TRUE) {
  abs(x - mean(x, na.rm = na.rm)) <= thres * sd(x, na.rm = na.rm)
}
isnt_out_mad <- function(x, thres = 3, na.rm = TRUE) {
      abs(x - median(x, na.rm = na.rm)) <= thres * mad(x, na.rm = na.rm)
}
isnt_out_tukey <- function(x, k = 1.5, na.rm = TRUE) {
  quar <- quantile(x, probs = c(0.25, 0.75), na.rm = na.rm)
  iqr <- diff(quar)
  
  (quar[1] - k * iqr <= x) & (x <= quar[2] + k * iqr)
}
maha_dist <- . %>% select_if(is.numeric) %>%
    mahalanobis(center = colMeans(.), cov = cov(.))

isnt_out_maha <- function(tbl, isnt_out_f, ...) {
  tbl %>% maha_dist() %>% isnt_out_f(...)
}
isnt_out_funs <- funs(
  z = isnt_out_z,
  mad = isnt_out_mad,
  tukey = isnt_out_tukey
)

## Warning: funs() is soft deprecated as of dplyr 0.8.0
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once per session.

dt1_imp <- dt1 %>% 
  drop_na(cluster_vars) %>%
  cbind(rock = as.factor(kmeans(dt1_knn, centers = 8)$cluster)) %>% 
  filter(!(abs(F80_M1 - median(F80_M1)) > 4*sd(F80_M1))) %>%
  filter(!(abs(hard - median(hard)) > 5*sd(hard))) %>% 
  filter(FF_1 > 0)
  
# intervals <- dt1_imp %>%
#   group_by(rock) %>%
#   summarise(min = min(hard), max = max(hard))
# intervals

ggpairs(dt1_imp %>% select(cluster_vars, rock),
      columns= 1:(1 + length(cluster_vars)), aes(colour = rock, alpha=0.5),
      lower=list(continuous="points"),
      upper=list(continuous="blank"),
      axisLabels="none", switch="both") +
      theme_bw()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

There is some potential with bivariate clustering. A more statisticall and technically correct approach would be to generate samples and test hypothesis with ANOVA.

5.2 Breaking harness into categories

ApplyTerciles <- function(x) {
  cut(x, breaks = c(quantile(dt1_imp$hard, 
                             probs = seq(0, 1, by = 1/3),
                             na.rm = TRUE)), 
      labels=c("soft","medium","hard"), 
      include.lowest=TRUE)}

dt1_imp$hard_qt <- sapply(dt1_imp$hard, ApplyTerciles)

ApplyQuintiles <- function(x) {
  cut(x, breaks = c(quantile(dt1_imp$F80_M1, 
                           probs = seq(0, 1, by = 0.2),
                           na.rm = TRUE)), 
  # cut(x, breaks = c(80,90,100,110, 120,Inf), 
        labels=c("fine_1","fine_2","fine_3", 
                 "coarse_1", "coarse_2"), 
      include.lowest=TRUE)}
dt1_imp$F80_qt <- sapply(dt1_imp$F80_M1, ApplyQuintiles)


p1 <- dt1_imp %>%
  mutate(visitors = hard, wday = hard_qt) %>% 
  group_by(wday) %>%
  summarise(mean_log_visitors = mean((visitors)),
            sd_log_visitors = sd((visitors))) %>%
  ggplot(aes(wday, mean_log_visitors, color = wday)) +

  geom_point(size = 4) +
  geom_errorbar(aes(ymin = mean_log_visitors - sd_log_visitors,
                    ymax = mean_log_visitors + sd_log_visitors,
                    color = wday), width = 0.5, size = 0.7) +
  theme(legend.position = "none")

p2 <- dt1_imp %>%
  mutate(visitors = hard, wday = hard_qt) %>% 
  mutate(visitors = log1p(visitors)) %>%
  ggplot(aes(visitors, wday, fill = wday)) +
  geom_density_ridges(bandwidth = 0.1) +
  theme(legend.position = "none") +
  labs(x = "log1p(visitors)", y = "")

p3 <- dt1_imp %>%
  mutate(visitors = F80_M1, month = F80_qt) %>% 
  group_by(month) %>%
  summarise(mean_log_visitors = mean(log1p(visitors)),
            sd_log_visitors = sd(log1p(visitors))) %>%
  ggplot(aes(month, mean_log_visitors, color = month)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = mean_log_visitors - sd_log_visitors,
                    ymax = mean_log_visitors + sd_log_visitors,
                    color = month), width = 0.5, size = 0.7) +
  theme(legend.position = "none")

p4 <- dt1_imp %>%
  mutate(visitors = F80_M1, month = F80_qt) %>% 
  mutate(visitors = log1p(visitors)) %>%
  ggplot(aes(visitors, month, fill = month)) +
  geom_density_ridges(bandwidth = 0.1) +
  theme(legend.position = "none") +
  labs(x = "log1p(visitors)", y = "")


grid.arrange(p1, p2, p3, p4, nrow=2, ncol=2)

5.3 Effects of hardness

We try to see if the hardness has any effects on the efficiency defined as tph and operating bond wond index.

petal.w.l <- ggplot(data = dt1_imp, mapping = aes(x = W_1, y = F_M1, color = F80_qt)) +
      geom_point() +
      geom_smooth() +
      scale_color_brewer(palette = 'Accent') +
      theme_classic() +
      xlim(4.5, 10) +
      ylim(1000, 4000) +
      theme(plot.background = element_rect(fill = "grey97")) +
      labs(title = 'Effects of Feed size', 
           x = 'Bond Operating Index', 
           y = 'Mill Throughput (TpH)')
ggplotly(petal.w.l)

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 2416 rows containing non-finite values (stat_smooth).

No evident correlation.

6 Start up Effect

6.1 Anomaly detection with Isolation Forest

import pandas as pd
import numpy as np
from sklearn.ensemble import IsolationForest
from sklearn.preprocessing import StandardScaler
import matplotlib
import matplotlib.pyplot as plt
plt.style.use('ggplot')
  
df1_imp = r.dt1_imp.copy()
df1_imp = df1_imp.assign(
    date_time = pd.to_datetime(df1_imp.loc[:, "date"])
    ).loc[:, ["date_time", "V_M1"]].sort_values('date_time')

df1_imp = df1_imp.assign(
  date_time_int = df1_imp.loc[:, "date_time"].astype(np.int64)
  ).loc[:, ["date_time_int", "V_M1"]]
# df1_imp

outliers_fraction = 0.05
scaler = StandardScaler()
np_scaled = scaler.fit_transform(df1_imp)
data = pd.DataFrame(np_scaled)
# train isolation forest
model =  IsolationForest(behaviour='new', contamination = outliers_fraction)
model.fit(data)

## IsolationForest(behaviour='new', bootstrap=False, contamination=0.05,
##                 max_features=1.0, max_samples='auto', n_estimators=100,
##                 n_jobs=None, random_state=None, verbose=0, warm_start=False)

df1_imp = df1_imp.assign(
anomaly = pd.Series(model.predict(data)))

fig, ax = plt.subplots(figsize=(10,6))
a = df1_imp.loc[df1_imp["anomaly"] == - 1, ['date_time_int', 'V_M1']]
ax.plot(df1_imp['date_time_int'], df1_imp['V_M1'], color='blue', label = 'Normal')
ax.scatter(a['date_time_int'], a['V_M1'], color = 'red', label = 'Anomaly')
plt.legend()
plt.show()

6.2 Velocity trend over time

dts <- dt1_imp %>% 
  cbind(py$df1_imp %>% select(c(anomaly))) %>% 
  mutate(anomaly = ifelse(anomaly == -1, 0, anomaly))

dts$Date <- as.Date(dts$date, "%Y-%m-%d")

dts <- dts[order(as.Date(dts$Date, format="%Y-%m-%d")),]

dts %>%
  group_by(Date) %>%
  summarise(anomaly = sum(anomaly)) %>%
  ggplot(aes(Date, anomaly)) +
  geom_line() +
  geom_smooth(method = "loess", color = "blue", span = 1/7) +
  labs(y = "Velocity SAG 1 (RPM)", x = "Date")

Overall, we can see a decreasing trend in Velocity anomalies.

7 Fresh Feed Trend

7.1 Time series data

# dts$Date <- as.Date(dts$date, "%Y-%m-%d")

# dts <- dts[order(as.Date(dts$Date, format="%Y-%m-%d")),]

dts <- dt1_imp %>% 
  mutate(Date = as.Date(dt1_imp$date, "%Y-%m-%d"),
         Date = ymd(Date),
         year = as.factor(year(Date)),
         color = ifelse(FF_1 > mean(.$FF_1), 'Above', 'Below'),
         month = month(Date),
         day = day(Date)) %>% 
  cbind(py$df1_imp %>% select(c(anomaly))) %>% 
  # filter(anomaly == 1) %>% 
  select(Date, year, month, day, 
         color, FF_1, F80_qt, hard_qt,
         F_M1, Ecs_1, CC_M1)
headTail(dts)

##             Date year month day color    FF_1   F80_qt hard_qt    F_M1
## 0     2018-01-01 2018     1   1 Above 2496.15   fine_2    hard 2947.63
## 1     2018-01-01 2018     1   1 Above 2457.71   fine_1    hard 2892.32
## 2     2018-01-01 2018     1   1 Below 2104.06   fine_3    hard 2605.67
## 3     2018-01-01 2018     1   1 Above 2664.65   fine_3    hard 2840.89
## ...         <NA> <NA>   ... ...  <NA>     ...     <NA>    <NA>     ...
## 13378 2019-08-26 2019     8  26 Below 1996.02 coarse_1  medium 1996.02
## 13379 2019-08-26 2019     8  26 Below 1974.99   fine_3  medium 1974.99
## 13380 2019-08-26 2019     8  26 Below 1501.83   fine_2  medium 1501.83
## 13381 2019-08-26 2019     8  26 Below 1826.39   fine_3  medium 1826.39
##       Ecs_1   CC_M1
## 0      4.83 4199.25
## 1      4.99 4197.64
## 2      5.52 4189.02
## 3      5.06 4190.68
## ...     ...     ...
## 13378  6.13 3523.25
## 13379  5.99 3528.32
## 13380  7.13 3472.14
## 13381  6.64 3572.53

7.2 Fresh Mill Throughput yearly trend

ggplot(data = dts, aes(x = year, y = FF_1)) +
        geom_jitter(pch = 21, alpha = .2, color = 'dark orange') +
        geom_boxplot(color = 'dark blue') +
        theme_few() +
        theme(legend.position = 'none') +
        geom_hline(yintercept = mean(dts$FF_1), linetype = 'dotted') +
        xlab('') +
        ylab('Mill Throughput (TpH)') +
        labs(title = 'Yearly performance',
             subtitle = 'Mill Throughput in 2018-2019')

7.3 Fresh Mill Throughput monthly trend

ggplot(data = dts, aes(x = as.factor(month), y = FF_1)) +
        geom_jitter(pch = 21, alpha = .2, color = 'dark orange') +
        geom_boxplot(color = 'dark blue') +
        theme_few() +
        theme(legend.position = 'none') +
        geom_hline(yintercept = mean(dts$FF_1), linetype = 'dotted') +
        xlab('Month') +
        ylab('Mill Throughput (TpH)') +
        labs(title = 'Monthly performance',
             subtitle = 'Mill Throughput in 2018-2019')

7.4 Fresh Mill Throughput trend

dts %>%
  group_by(Date) %>%
  summarise(FF_1 = mean(FF_1)) %>%
  ggplot(aes(Date, FF_1)) +
  geom_line() +
  geom_smooth(method = "loess", color = "blue", span = 1/7) +
  labs(y = "Mill Thoughput (TpH)", x = "Date")

7.5 Fresh Mill Throughput yearly trend by feed size

ggplot(data = dts, aes(x = Date, y = FF_1, color = F80_qt)) + 
  geom_line(size = 0.5) +
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC"),
        legend.text = element_text(size=7)) + 
  labs(title="Fresh Mill Throughput by feed size") + 
  geom_hline(yintercept=max(dts$FF_1), linetype="dashed", color = "red") +
  geom_hline(yintercept=min(dts$FF_1), linetype="dashed", color = "blue")

7.6 Fresh Mill Throughput yearly trend by hardness

ggplot(data = dts, aes(x = Date, y = FF_1, color = hard_qt)) + 
  geom_line(size = 0.5) +
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC"),
        legend.text = element_text(size=7)) + 
  labs(title="Fresh Mill Throughput by hardness") + 
  geom_hline(yintercept=max(dts$FF_1), linetype="dashed", color = "red") +
  geom_hline(yintercept=min(dts$FF_1), linetype="dashed", color = "blue")

7.7 Fresh Mill Throughput yearly CUMSUM approach

F80_lst = as.list(unique(as.character(dts$year)))

for (i in seq_along(F80_lst)) {
  summer_center <- mean(dts$FF_1)
  summer_sd <- sd(dts$FF_1)
  F80_index <- F80_lst[[i]]
  df <- dts %>%
    filter(month < 8) %>% 
    filter(as.character(year) == F80_index) %>%
    dplyr::select(FF_1)
  # fit a cusum model to that year
  qsum <- qcc::cusum(
          data = df$FF_1,
          centervalue = summer_mean,
          std.dev = summer_sd,
          se.shift = 1.96,
          plot = T)
  }

7.8 Filtering by type

rock_1 <- dts %>% 
  filter(hard_qt == "soft")
rock_1 <- as_tbl_time(rock_1, index=Date)
rock_1 <- as_period(rock_1, '1 month')

rock_2 <- dts %>% 
  filter(hard_qt == "medium")
rock_2 <- as_tbl_time(rock_2, index=Date)
rock_2 <- as_period(rock_2, '1 month')

rock_3 <- dts %>% 
  filter(hard_qt == "hard")
rock_3 <- as_tbl_time(rock_3, index=Date)
rock_3 <- as_period(rock_3, '1 month')

7.9 Fresh Mill Throughput decomposition

# options(repr.plot.width=8, repr.plot.height=6)
rock_1_monthly <- rock_1 %>%
  ggplot(aes(x = Date, y = FF_1, color = F80_qt)) + 
  geom_line(size = 1.5) +
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC"),
        legend.text = element_text(size=7)) + 
  labs(title="Soft rock") +
  geom_hline(yintercept=max(rock_1$FF_1), linetype="dashed", color = "red") +
  geom_hline(yintercept=min(rock_1$FF_1), linetype="dashed", color = "blue")

# Let's create a volume chart
rock_1_volume <- rock_1 %>%
  ggplot(aes(x=Date, y = FF_1)) + 
  geom_bar(stat='identity', fill="#7FB3D5",
           color="black") + theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC")) +
  geom_smooth(method="loess", color="red")

rock_2_monthly <- rock_2 %>% 
  ggplot(aes(x=Date, y = FF_1, color = F80_qt)) + 
  geom_line(size = 1.5) +
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC"),
        legend.text = element_text(size=7)) + 
  labs(title="Medium") + 
  geom_hline(yintercept=max(rock_2$FF_1), 
             linetype="dashed", color = "red") + 
  geom_hline(yintercept=min(rock_2$FF_1), linetype="dashed", color = "blue")

rock_2_volume <- rock_2 %>%
  ggplot(aes(x=Date, y=FF_1)) +
  geom_bar(stat='identity',fill="#58D68D",color="black") + 
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC")) +
  geom_smooth(method="loess", color="red")

rock_3_monthly <- rock_3 %>% 
  ggplot(aes(x=Date, y=FF_1, color=F80_qt)) + 
  geom_line(size = 1.5) +
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC"),
        legend.text = element_text(size=7)) + 
  labs(title="Hard rock") + 
  geom_hline(yintercept=max(rock_3$FF_1), 
             linetype="dashed", color = "red") + 
  geom_hline(yintercept=min(rock_3$FF_1), linetype="dashed", color = "blue")

rock_3_volume <- rock_3 %>%
  ggplot(aes(x=Date, y=FF_1)) +
  geom_bar(stat='identity',fill="#FE642E",color="black") + 
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC")) +
  geom_smooth(method="loess", color="red")

grid.arrange(rock_1_monthly, rock_2_monthly, rock_3_monthly,
             rock_1_volume, rock_2_volume, rock_3_volume, 
             nrow=2, ncol=3)

7.10 Load Cell

7.10.1 Load trend

dts %>%
  group_by(Date) %>%
  summarise(CC_M1 = mean(CC_M1)) %>%
  ggplot(aes(Date, CC_M1)) +
  geom_line() +
  geom_smooth(method = "loess", color = "blue", span = 1/7) +
  labs(y = "Mill Thoughput (t)", x = "Date")

7.10.2 Load yearly trend by feed size

ggplot(data = dts, aes(x = Date, y = CC_M1, color = F80_qt)) + 
  geom_line(size = 0.5) +
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC"),
        legend.text = element_text(size=7)) + 
  labs(title="Load cell by feed size") + 
  geom_hline(yintercept=max(dts$CC_M1), linetype="dashed", color = "red") +
  geom_hline(yintercept=min(dts$CC_M1), linetype="dashed", color = "blue")

7.10.3 Load yearly trend by hardness

ggplot(data = dts, aes(x = Date, y = CC_M1, color = hard_qt)) + 
  geom_line(size = 0.5) +
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC"),
        legend.text = element_text(size=7)) + 
  labs(title="Load Cell by hardness") + 
  geom_hline(yintercept=max(dts$CC_M1), linetype="dashed", color = "red") +
  geom_hline(yintercept=min(dts$CC_M1), linetype="dashed", color = "blue")

7.10.4 Load yearly CUMSUM approach

F80_lst = as.list(unique(as.character(dts$year)))

for (i in seq_along(F80_lst)) {
  summer_center <- mean(dts$CC_M1)
  summer_sd <- sd(dts$CC_M1)
  F80_index <- F80_lst[[i]]
  df <- dts %>%
    filter(month < 8) %>% 
    filter(as.character(year) == F80_index) %>%
    dplyr::select(CC_M1)
  # fit a cusum model to that year
  qsum <- qcc::cusum(
          data = df$CC_M1,
          centervalue = summer_mean,
          std.dev = summer_sd,
          se.shift = 1.96,
          plot = T)
  }

7.10.5 Load Cell decomposition

# options(repr.plot.width=8, repr.plot.height=6)
rock_1_monthly <- rock_1 %>%
  ggplot(aes(x = Date, y = CC_M1, color = F80_qt)) + 
  geom_line(size = 1.5) +
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC"),
        legend.text = element_text(size=7)) + 
  labs(title="Soft rock") +
  geom_hline(yintercept=max(rock_1$CC_M1), linetype="dashed", color = "red") +
  geom_hline(yintercept=min(rock_1$CC_M1), linetype="dashed", color = "blue")

# Let's create a volume chart
rock_1_volume <- rock_1 %>%
  ggplot(aes(x=Date, y = CC_M1)) + 
  geom_bar(stat='identity', fill="#7FB3D5",
           color="black") + theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC")) +
  geom_smooth(method="loess", color="red")

rock_2_monthly <- rock_2 %>% 
  ggplot(aes(x=Date, y = CC_M1, color = F80_qt)) + 
  geom_line(size = 1.5) +
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC"),
        legend.text = element_text(size=7)) + 
  labs(title="Medium") + 
  geom_hline(yintercept=max(rock_2$CC_M1), 
             linetype="dashed", color = "red") + 
  geom_hline(yintercept=min(rock_2$CC_M1), linetype="dashed", color = "blue")

rock_2_volume <- rock_2 %>%
  ggplot(aes(x=Date, y=CC_M1)) +
  geom_bar(stat='identity',fill="#58D68D",color="black") + 
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC")) +
  geom_smooth(method="loess", color="red")

rock_3_monthly <- rock_3 %>% 
  ggplot(aes(x=Date, y=CC_M1, color=F80_qt)) + 
  geom_line(size = 1.5) +
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC"),
        legend.text = element_text(size=7)) + 
  labs(title="Hard rock") + 
  geom_hline(yintercept=max(rock_3$CC_M1), 
             linetype="dashed", color = "red") + 
  geom_hline(yintercept=min(rock_3$CC_M1), linetype="dashed", color = "blue")

rock_3_volume <- rock_3 %>%
  ggplot(aes(x=Date, y=CC_M1)) +
  geom_bar(stat='identity',fill="#FE642E",color="black") + 
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5),
        plot.background=element_rect(fill="#D5D8DC")) +
  geom_smooth(method="loess", color="red")

grid.arrange(rock_1_monthly, rock_2_monthly, rock_3_monthly,
             rock_1_volume, rock_2_volume, rock_3_volume, 
             nrow=2, ncol=3)